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A computer program is presented by which one may calculate the mul- 
tiple electric dipole, electric quadrupole and magnetic dipole Coulomb ex- 
citation with relativistic heavy ions. The program applies to an arbitrary 
nucleus, specified by the spins and energies of the levels and by all El, E2 
and Ml matrix elements. Nuclear excitation is calculated optionally for 
monopole, dipole and quadrupole excitations and needs inputs of optical 
potentials. For given bombarding conditions, the differential cross sec- 
tions and statistical tensors (useful to calculate 7-ray angular distribution 
functions) are computed. 

1 Introduction 

Relativistic Coulomb Excitation (RCE) is a well established tool to unravel inter- 
esting aspects of nuclear structure [l|, Q]. The RCE induced by large-Z projectiles 
and/or targets, often yields large cross sections in grazing collisions. This results 
from the large nuclear response (specially in the region of the giant resonances) to 
the acting electromagnetic fields. As a consequence, a strong coupling between the 
excited states is expected. The present report describes a computer program for 
the calculation of multiple excitation among a finite number of nuclear states. The 
system of coupled differential equations for the time-dependent amplitudes of the 
eigenstates of the free nucleus is solved numerically for electric dipole (El), elec- 
tric quadrupole (E2), and magnetic dipole (Ml) excitations. If the optical potential 
for the system is given, the program can also calculate the amplitudes for (nuclear) 
monopole, dipole and quadrupole excitations. All nuclear quantities, either known 
from experiments or calculated from a model, as well as the conditions realized in 
the experiment, are explicitly specified as input parameters. The program then com- 
putes the Coulomb + nuclear excitation probabilities and cross sections as well as the 
statistical tensors for the angular distribution of the 7-quanta. 



2 The semiclassical method for the CC-problem 

In relativistic heavy ion collisions, the wavelength associated to the projectile-target 
relative motion is much smaller than the characteristic lengths of the system. It is, 
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therefore, a reasonable approximation to treat r as a classical variable r(t), given 
at each instant by the trajectory followed by the relative motion. At high energies 
it is also a good approximation to replace this trajectory by a straight line. The 
intrinsic dynamics can then be handled as a quantum mechanics problem with a 
time dependent Hamiltonian. This treatment is discussed in full details by Alder and 
Winther in ref. 0. 

The intrinsic state \tp(t) > satisfies the Schrodinger equation 

{H + V [rm\m) = i*^P-- (!) 

Above, H Q is the intrinsic Hamiltonian and V is the channel-coupling interaction. 

Expanding the wave function in the set {| k); k = 1,N} of eigenstates of H , 
where N is the number states included in the coupled-channels (CC) problem, we 
obtain a set of coupled equations. Taking the scalar product with each of the states 
< j\, we get 

N 

ih a k (t) = Y, (k\ V(t) | j) exp [i(E k - EJt/h] a^t) , k, j = 1 to N . (2) 

where E n is the energy of the state |n) . It should be remarked that the amplitudes 
depend also on the impact parameter b specifying the classical trajectory followed 
by the system. For the sake of keeping the notation simple, we do not indicate 
this dependence explicitly. We write, therefore, a n (t) instead of a n (b,t). Since the 
interaction V vanishes as t — > ±00, the amplitudes have as initial condition a n (t — > 
—00) = 5 n i and they tend to constant values as t — > 00. 

A convenient measure of time is given by the dimensionless quantity r = jvt/b, 
where 7 = (1 — v 2 /c 2 ) -1 / 2 is the Lorentz factor for the projectile velocity v. A 
convenient measure of energy is E Q = 'jhv/b. In terms of these quantities the CC 
equations become 



da dr) .^, II „ MI , ,. t \ , , x V{t) t E k — E 



j 



di~ j= i &o ^0 

(3) 

The nuclear states are specified by the spin quantum numbers I and M. Therefore, 
the excitation probability of an intrinsic state | n) =| J n , M n ) in a collision with impact 
parameter b is obtained from an average over the initial orientation (Mi) , and a sum 
over the final orientation of the nucleus, respectively: 

p n(b) = ^—r E K^MOI 2 . (4) 

ZI 1 + 1 Mi,M n 

The total cross section for excitation of the state \n > is obtained by the classical 
expression 

a n = 2tt J P n {b) b db . (5) 
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2.1 Coulomb excitation 



We consider a nucleus 2 which is at rest and a relativistic nucleus 1 which moves along 
the z-axis. Nucleus 2 is excited from the initial state \IjMj > to the state \IkMk > 
by the electromagnetic field of nucleus 1. The nuclear states are specified by the 
spin quantum numbers Ij, Ik and by the corresponding magnetic quantum numbers 
Mj and We assume that the relativistic nucleus 1 moves along a straight-line 
trajectory with impact parameter b, which is therefore also the distance of the closest 
approach between the center of mass of the two nuclei at the time t = 0. The 
interaction, Vc(t), due to the electromagnetic field of the nucleus 1 acting on the 
charges and currents nucleus 2 can be expanded into multipoles, as explained in ref. 
|§. One has 

^c(r) = ^ = £WWr) , (6) 



E() 



where n = E, M denotes electric and magnetic interactions, respectively, and (mis- 
prints in ref. 0] have been corrected) 



W nXfl (r 



i Z\£ i 



2tt 



v ' hvb x \ V (2A + 1)!! 
where Ai(n\,fi) is the multipole moment of order A/i 



Qir\»{t,r) M{ir\,ii) 



(7) 



M(E\,im)= d z r p{v) r A F lM (r) , 



and 



M(Ml,fi) = ~y c J dh J(r)-L(ry lM ) 



(8) 



(9) 



p (J) being the nuclear charge (current). The quantities Q-kx^t) were calculated in 
ref. H, and for the El, E2, and Ml multipolarities are given by 



<3i?io(£, t) = 7 V2 



r0 3 (r)-^Q 2 0(r) 



Qei±i(£,t) 



Qmvo(£, r) = ; Q M i±i{£, t) = i[ - ) <p :i (r) , 



and 



(10) 
(11) 



Qe2o(£,t) = 7 \ 
Qe2±i(€,t) = ±7 



(2r 2 -l) 5 (r)-^ (l) 2 r<f> 3 (r) 
6r0 5 (r)-^(^) 2 3 (r) ; Q S2±2 (r) 



(r), (12) 



3 



where <fi (r) = (1 + r 2 )~ / . 

We use here the notation of Edmonds |J where the reduced multipole matrix 
element is defined by 

M kj {E\n) = (-1)'*-"* ( J^*^ J < hWMiEX)^ > (13) 

To simplify the expression (§) we introduce the dimensionless parameter ip^ by the 
relation 

Then we may write eq. |3] in the form 

dak j- T ^ =-iJ2 J2Q^^kj,r)ip ( k f exp(i6yr) a^r) . (14) 

The fields QttXh(^ t ) peak around r = 0, and decrease rapidly within an interval 
At ~ 1, corresponding to a collisional time At ~ b/yv. This means that numerically 
one needs to integrate the CC equations in time within an interval of range n x At 
around r = 0, with n equal to a small integer number. 



2.2 Nuclear excitation 

In peripheral collisions the nuclear interaction between the ions can also induce excita- 
tions. According to the Bohr-Mottelson particle-vibrator coupling model, the matrix 
element for the transition j — > k is given by 

=< >= < IMY^IjMj > Y Xfi (v)U x (r) 

(15) 

where S x is the vibrational amplitude and U\(r) is the transition potential. 

The transition potentials for nuclear excitations can be related to the optical 
potential in the elastic channel. This is discussed in details in ref. |J. The transition 
potentials for isoscalar excitations are 

U (r) = 3U opt (r) + , (16) 



for monopole 
for dipole, and 



tfi(r) = ^ + i Ro^ft , (it) 
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(r) = dU opt (r) 



dr 

for quadrupole modes. 

The time dependence of the matrix elements above can be obtained by making a 
Lorentz boost. One gets 

V$$L(t) = <I k M k \U\I j M j > 



5> 
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V2A + 1 



< IkMklYx^Mj > Y Xfl (B(t), 0) U x [r{t)} , (19) 



where r(t) = \/b 2 + / y 2 v 2 t 2 = 1/ (b(j) (r)) , 9 = t0(t), and 



< hM^Yx^Mj >= (-1 



Jk~M k 



(24 + l)(2A+l) 
4vr(2i j + 1) 



1/2 



Ik A Ij \ I Ik X Ij 
-MkuMj) \ 



(20) 

To put it in the same notation as in eq. fll^D , we define Q^\\{t) = V^ x ^(t) / E Q , 
and the coupled-channels equations become 



da k (j) 
dr 



N 

-EE 

j=l A/it 



Q 



exp(i£ kj r) aj(r) 



(21) 



2.3 Absorption at small impact parameters 

If the optical potential U opt (r) is known, the absorption probability in grazing colli- 
sions can be calculated in the eikonal approximation as 



A(b) = exp 



2 r°° 

— I lm[U opt (r)]dz 



(22) 



where r = \/b 2 + z 2 . If the optical potential is not known, the absorption probability 
can be calculated from the optical limit of the Glauber theory of multiple scattering, 
which yields: 



A(b) = exp j-(TArjv J J Pi(r') p 2 (r - r) d 3 r' 



dz 



(23) 



where <Jnn = 40 mb is the nucleon-nucleon cross section in high energy collisions 
and pi is the ground state density of the nucleus i. These densities are taken from 
the droplet model densities of Myers and Swiatecki || , but can be easily replaced by 
more realistic densities. 

Including absorption, the total cross section for excitation of the state \n > is 
obtained by 

o n = 2ir J A{b)P n {b) bdb . (24) 
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2.4 Coulomb recoil correction 

At intermediate bombarding energies a correction of the excitation amplitudes due to 
Coulomb recoil can be easily introduced. The excitation occurs most probably when 
the nuclei are closer along the trajectory. But then, due to recoil, they are displaced 
by an extra distance of order of a /7, where a = Z 1 Z 2 e 2 /m v 2 . It is thus expected 
that one should correct the amplitudes by a rescaling of the impact parameter. In 
ref. U] it was shown that a comparison between the Rutherford integrals with large 
angular momenta for the non-relativistic case with the straight-line integrals for the 
relativistic case yields a reasonable correction. It amounts in the replacement of b of 
P n {b) in equation fl2~j|) (and also in A{b) of eq. (|2"2"D) by b' = b + TTa /2j, i.e., for a 
given impact parameter b the excitation amplitudes and probabilities are calculated 
for b'. 

2.5 Gamma-ray angular distributions 

As for the non-relativistic case || [3|, the angular distributions of gamma rays following 
the excitation depend on the frame of reference used. In our notation, the z-axis 
corresponds to the beam axis, and the statistical tensors are given by 



AUN) ~ Wi + 1) ^ , \M N M' n k 

V 1 ^ M n =-(M' n +k),M n V N ' 

x E fl UW fl W^). (25) 

Mi 

where N is the state from which the gamma ray is emitted, and denotes the initial 
state of the nucleus, before the excitation. To calculate the angular distributions of 
the gamma rays one needs the statistical tensors for k = 0, 2, 4 and — k < k < k. Since 
the nuclear levels are not only populated by Coulomb excitation but also by internal 
conversion and gamma transitions cascading down from higher states, the angular 
distributions depend not only on the 7-transition rates, S^^m (f° r ah M's), but 
also on the internal conversion coefficients, a^-,M- With those parameters at hand 
one calculates the angular distributions of gamma rays, dW^^^M / dVL^, by using the 
equations given in section IV of ref. (also reprinted in [[|). 



3 Computer program and user's manual 

The units used in the program are fm (Fermis) for distances and MeV for energies. 
The output cross sections are given in millibarns. 
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3.1 Input parameters 

To avoid exceeding use of computer's memory, the file RELEX.DIM contains the di- 
mension of the arrays and sets in the maximum number of levels (NMAX) , maximum 
total number of magnetic substates, (NSTMAX), maximum number of impact pa- 
rameters (NBMAX), and maximum number of coordinates points used in the optical 
potentials and absorption factors, (NGRID). 

The integrals appearing in eqs. (22-24) are performed by the 1/3-Simpson's in- 
tegration rule. It is required that NGRID be a even number, since an extra point 
(origin) is generated in the program. 

The file RELEX.IN contains all other input parameters. These are 

1. AP, ZP, AT, ZT, which are the projectile and the target mass and charge num- 
bers, respectively. 

2. ECA, the bombarding energy per nucleon in MeV. 

3. EX(j) and SPIN(j): the energy and spins of the individual states j. 

4. MATEl(j,k), MATE2(j,k), MATMl(j,k), the reduced matrix elements for El, 
E2 and Ml excitations, j k, (as defined in (|i~3|)). 

5. DELTE(OJ), DELTE(lj), DELTE(2,j), the deformation parameters for monopole, 
dipole and quadrupole nuclear excitations, entering eq. (|19|) . 

The input cards in file RELEX.IN are organized as following: 

1. First data card: AP, ZP, AT, ZT, ECA, IW, IOUT. 
IW = (or 1) for projectile (or target) excitation 

IOUT = (or 1) for output (or not) of statistical tensors. The statistical 
tensors are calculated for each impact parameter, so that one can use them in 
the computation of o?P 7 7v^m(^) / 'dQ, y = P^{b) .dW^N-*M I f d£l 7 . 

2. Second data card: NB, ACCUR, BMIN, ITOT 
NB = number of points in impact parameter mesh. 

ACCUR = accuracy required for the time integration of the CC-equations for 
each impact parameter. A reasonable value is ACCUR = 0.001, i.e., 0.1%. 

BMIN = minimum impact parameter (can be set to zero). 

ITOT = 1 (0) prints (does not print) out impact parameter probabilities. 

3. Third data card: IOPW, IOPNUC 

IOPW is an option to use (or not) an optical potential: IOPW = 1 (0). 
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If the optical potential is provided (IOPW = 1), it should be stored in file 
IOPW.IN in rows of R x Real[U(R)] x Imag[U(R)]. The first row in this file gives 
the number of rows (maximum = NGRID). The program makes an interpolation 
to obtain intermediate values. 

IOPNUC = 1 (0) is an option to compute (or not) nuclear excitation. 

4. Fourth data card: NST 
Number of nuclear levels. 

5. Fifth and following data cards: J, EX(J), SPIN(J) 

Input of state labels (J), energy (EX), and angular momentum (SPIN). J ranges 
from 1 to NST and should be listed in increasing value of energies. 

6. Following data cards: J, K, MATE1(J,K), MATE2(J,K), MATM1(J,K) 

Reduced matrix elements for electromagnetic El, E2 and Ml excitations: 

MATEl(j^k), for electric dipole transitions, MATE2(j^k), for electric quadrupole 
transitions, MATMl(j— >k), for magnetic dipole transitions. Matrices for reori- 
entation effects, i— >i, can also be given. 

To end reading the matrix elements add a row of zeros for each column. 

7. Following data cards: K, DELTE (0,K), DELTE(1,K), DELTE(2,K) 

If IOPNUC = 1 these rows give the values of deformation parameters for 
monopole, dipole and quadrupole nuclear excitations, respectively, for each ex- 
cited state K. 

3.2 Computer program 

The program starts with a catalogue of the nuclear levels by doing a correspondence 
of integers to each magnetic substate. J = 1 corresponds to the lowest energy level, 
with the magnetic quantum number M 1 = J increases with M 1 and so on for 

the subsequent levels. 

A mesh in impact parameter is done, reserving half of the impact parameter 
points, i.e., NB/2, to a finer mesh around the grazing impact parameter, defined as 
bo = 1.2 (^A 1 / 3 + A^ 3 ^J fm. The interval bo/2 fm < b < 3bo/2 fm is covered by this 
mesh. A second mesh, with the other half of points, extends from b = 3bo/2 fm to 
b = 200 fm. Except for the very low excitation energies (Ej <C 1 MeV), combined 
with very large bombarding energies (7 ^> 1), this upper value of b corresponds to 
very small excitation probabilities, and the calculation can be safely stopped. The 
reason for a finer mesh at small impact parameters is to get a good accuracy at the 
region where both nuclear, Coulomb, and absorption factors play equally important 
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roles. At large impact parameters the probabilities fall off smoothly with b, justifying 
a wider integration step. 

A mesh of NGRID points in polar coordinates is implemented to calculate the 
nuclear excitation potentials and absorption factors, according to the equations pre- 
sented in sections 2.2 and 2.3. The first and second derivatives of the optical potentials 
are calculated by the routine DERIVATIVE. A 6-point formula is used for the pur- 
pose. The routines TWOFOLD computes the folding over the densities, as used in 
eq. (p3|) . Routines RHOPP and RHONP generate the liquid drop densities, and the 
routine PHNUC computes the eikonal integral appearing in eq. (|2"2"|). 

Repeated factors for the nuclear and for the Coulomb potentials are calculated 
in the main program and stored in the main program with the arrays PSITOT and 
PSINUC. These arrays are carried over in a common block to the routine VINT which 
computes the function Qn\\(t) + ^ Q^x^kj, T )^kj i use< ^ in e 1- (HI)- 

The time integrals are performed by means of an adaptive Runge-Kutta method. 
All routines used for this purpose have been taken from the Numerical Recipes Soft- 
ware, described in the book M. They are composed by the routines ODEINT, RKQS, 
RKCK, and RK4. The routine ODEINT varies the time step sizes to achieve the de- 
sired accuracy, controlled by the input parameter ACCUR. The right side of ( PH ) is 
computed in the routine DCADT, used externally by the fourth-order Runge-Kutta 
routine RK4. RKCK is a driver to increase time steps in RK4, and RKQS is used 
in ODEINT for the variation of step size and accuracy control. The main program 
returns a warning if the summed errors for all magnetic substates is larger than 10 x 
ACCUR. 

The routine THREEJ computes Wigner-3J coefficients (and Clebsh-Gordan coef- 
ficients), and YLM is used to compute the spherical harmonics. 

The routines SPLINE and SPLINT perform a spline interpolation of the excitation 
amplitudes, before they are used for integration by means of the routines QTRAP 
and QSIMP. All of them are from Numerical Recipes Software ||. 

The program delivers outputs in RELEX.OUT. The statistical tensors are given 
in RELEX2.0UT for each impact parameter (only if IOUT=l). 

The program is available at [http: / / www.if.ufrj .br^bertu/ prog. html . 
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